package info.openrocket.core.util;

import java.util.Random;

/**
 * A class that provides a source of pink noise with a power spectrum density
 * proportional to 1/f^alpha. The values are computed by applying an IIR filter
 * to
 * generated Gaussian random numbers. The number of poles used in the filter may
 * be
 * specified. Values as low as 3 produce good results, but using a larger number
 * of
 * poles allows lower frequencies to be amplified. Below the cutoff frequency
 * the
 * power spectrum density if constant.
 * <p>
 * The IIR filter use by this class is presented by N. Jeremy Kasdin,
 * Proceedings of
 * the IEEE, Vol. 83, No. 5, May 1995, p. 822.
 * 
 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
 */
public class PinkNoise {
	private final int poles;
	private final double[] multipliers;

	private final double[] values;
	private final Random rnd;

	/**
	 * Generate pink noise with alpha=1.0 using a five-pole IIR.
	 */
	public PinkNoise() {
		this(1.0, 5, new Random());
	}

	/**
	 * Generate a specific pink noise using a five-pole IIR.
	 * 
	 * @param alpha the exponent of the pink noise, 1/f^alpha.
	 */
	public PinkNoise(double alpha) {
		this(alpha, 5, new Random());
	}

	/**
	 * Generate pink noise specifying alpha and the number of poles. The larger the
	 * number of poles, the lower are the lowest frequency components that are
	 * amplified.
	 * 
	 * @param alpha the exponent of the pink noise, 1/f^alpha.
	 * @param poles the number of poles to use.
	 */
	public PinkNoise(double alpha, int poles) {
		this(alpha, poles, new Random());
	}

	/**
	 * Generate pink noise specifying alpha, the number of poles and the randomness
	 * source.
	 * 
	 * @param alpha  the exponent of the pink noise, 1/f^alpha.
	 * @param poles  the number of poles to use.
	 * @param random the randomness source.
	 */
	public PinkNoise(double alpha, int poles, Random random) {
		this.rnd = random;
		this.poles = poles;
		this.multipliers = new double[poles];
		this.values = new double[poles];

		double a = 1;
		for (int i = 0; i < poles; i++) {
			a = (i - alpha / 2) * a / (i + 1);
			multipliers[i] = a;
		}

		// Fill the history with random values
		for (int i = 0; i < 5 * poles; i++)
			this.nextValue();
	}

	public double nextValue() {
		double x = rnd.nextGaussian();
		// double x = rnd.nextDouble()-0.5;

		for (int i = 0; i < poles; i++) {
			x -= multipliers[i] * values[i];
		}
		System.arraycopy(values, 0, values, 1, values.length - 1);
		values[0] = x;

		return x;
	}

	// TODO: this method seems incomplete
	public static void main(String[] arg) {

		@SuppressWarnings("unused")
		PinkNoise source;

		source = new PinkNoise(1.0, 100);
		@SuppressWarnings("unused")
		double std = 0;
		for (int i = 0; i < 1000000; i++) {

		}

		// int n = 5000000;
		// double avgavg=0;
		// double avgstd = 0;
		// double[] val = new double[n];
		//
		// for (int j=0; j < 10; j++) {
		// double avg=0, std=0;
		// source = new PinkNoise(5.0/3.0, 2);
		//
		// for (int i=0; i < n; i++) {
		// val[i] = source.nextValue();
		// avg += val[i];
		// }
		// avg /= n;
		// for (int i=0; i < n; i++) {
		// std += (val[i]-avg)*(val[i]-avg);
		// }
		// std /= n;
		// std = Math.sqrt(std);
		//
		// System.out.println("avg:"+avg+" stddev:"+std);
		// avgavg += avg;
		// avgstd += std;
		// }
		// avgavg /= 10;
		// avgstd /= 10;
		// System.out.println("Average avg:"+avgavg+" std:"+avgstd);
		//
		// Two poles:

	}

}
